rm(list=ls())
library(foreign)
library(data.table)
library(RColorBrewer)
library(grid)
library(ggmap)
library(maptools)
library(dplyr)

options(scipen=10)

# shapefiles 

cz <- readShapePoly("CZ.shp")
cz.df <- fortify(cz)
cz.df <- data.frame(cz.df, cz@data[cz.df$id, ])
cz.df <- cz.df[cz.df$long < -50 & cz.df$long > -135 & cz.df$lat < 50 & cz.df$lat > 22, ]

names(cz.df)[8] <- "czone"

perc.rank <- function(x) trunc(rank(x))/length(x)

# reading czone-level data

cz_taa <- read.dta("kimpelc_czone.dta")
cz_taa$mean_l_TAAdlrs_wrks_cert_pc_qt <-  perc.rank(cz_taa$mean_l_TAAdlrs_wrks_cert_pc)
cz_taa$mean_l_trans_taaimp_pc_qt <-  perc.rank(cz_taa$mean_l_trans_taaimp_pc)
cz_taa$difference <- cz_taa$mean_l_TAAdlrs_wrks_cert_pc - cz_taa$mean_l_trans_taaimp_pc
cz_taa$difference_break <- cut(cz_taa$difference, breaks = c(-5505, -2000, -1500, -1000, -500, 0, 500, 1000, 1500, 2000, 17740),dig.lab=10)
cz_taa $success_rate <- cz_taa $certified_petitions_n/cz_taa $petitions_n

cz.df<- merge(cz.df, cz_taa, by = intersect(names(cz.df), names(cz_taa)))

# map base 

mapPlot <- ggplot(cz.df) + theme_bw() + theme(plot.background = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.border = element_blank() ,axis.title = element_blank(), axis.text = element_blank(), axis.ticks = element_blank(), legend.position="bottom",    legend.key.width = unit(1, "cm"), legend.text = element_text(size = 8))

# reading czone-by-decade data

cz.response <- read.dta("kimpelc_czone_t2.dta")
cz.response <- subset(cz.response, year!= 2007)

cz.response$exposure_qt <- perc.rank(cz.response$l_tradeusch_pw)
cz.response$response_qt <- perc.rank(cz.response$response)
cz.response$taa_qt <- perc.rank(cz.response$petitioningwrks_pw_mean)
cz.response$cert_response_qt <- perc.rank(cz.response$response_cert)
cz.response$response_qt9010 <- perc.rank(cz.response$response9010)

cz.response.1990 <- subset(cz.response, t2 == 0)
cz.response.2000 <- subset(cz.response, t2 == 1)

cz.change <- merge(cz.response.1990, cz.response.2000, by = c("czone"))

cz.change$exposure_delta <- cz.change$l_tradeusch_pw.y - cz.change$l_tradeusch_pw.x
cz.change$exposure_delta_qt <- perc.rank(cz.change$exposure_delta)

cz.change$response_delta <- cz.change$response.y - cz.change$response.x
cz.change$response_delta_qt <- perc.rank(cz.change$response_delta)

cz.change$taa_delta <- cz.change$petitioningwrks_pw_mean.y - cz.change$petitioningwrks_pw_mean.x
cz.change$taa_delta_qt <- perc.rank(cz.change$taa_delta)

cz.change$response_delta9010 <- cz.change$response9010.y - cz.change$response9010.x
cz.change$response_delta_qt9010 <- perc.rank(cz.change$response_delta9010)

cz.change <- cz.change[,c("czone", "response_delta", "exposure_delta", "taa_delta", "response_delta_qt", "exposure_delta_qt", "taa_delta_qt", "response_delta9010", "response_delta_qt9010")]

cz.df.1990 <- merge(cz.df, cz.response.1990, by = intersect(names(cz.df), names(cz.response.1990)))
cz.df.2000 <- merge(cz.df, cz.response.2000, by = intersect(names(cz.df), names(cz.response.2000)))

cz.df.change <- merge(cz.df, cz.change, by = intersect(names(cz.df), names(cz.change)))

# Appendix Figure A2 Difference in the estimated level of TAA payments based on TAA petitions versus UI

pdf("FigureA2.pdf", width=7, height=5)
mycolors = c(brewer.pal(n = 10, name = "RdBu"))

diff_map <- mapPlot + geom_polygon(data=cz.df, aes(x = long, y = lat, group = group, fill = cz.df$difference_break), lwd = 1/10, colour="lavenderblush4") + scale_fill_manual(name = "", values = mycolors) + theme(legend.text = element_text(size=8), legend.position="right", legend.key.size = unit(0.5, "cm"), legend.key.width=unit(0.5, "cm")) + guides(fill = guide_legend(ncol = 1, reverse = TRUE)) + theme(legend.margin=margin(l = -1.5, r = 1.5, unit='cm')) + coord_map(project="conic", lat0 = 30)

plot(diff_map)
dev.off()

# Appendix Figure A3 TAA Certification Rate, 1990-2007

pdf("FigureA3.pdf",width=7,height=5)

success_map <- mapPlot + geom_polygon(data=cz.df, aes(x = long, y = lat, group = group, fill = cz.df$success_rate), lwd = 1/10) +  scale_fill_continuous(name = "", guide="colorbar", limits=c(0.0,1.0)) +
  theme(legend.text = element_text(size=8),  legend.position="right" , legend.key.size = unit(1, "cm"), legend.key.width=unit(0.6, "cm")) + theme(legend.margin=margin(l = -1.5, r = 1.5, unit='cm')) + coord_map(project="conic", lat0 = 30)

plot(success_map)
dev.off()

# Appendix Figure A4 TAA Petitions and Certification across Commuting Zones 

pdf("FigureA4-A.pdf",width=5,height=5)

ggplot(cz_taa, aes(petitions_n, certified_petitions_n)) + xlab("Petitions, 1990-2007") + ylab("Certified Petitions, 1990-2007") + geom_point() + geom_smooth() + theme(text = element_text(size=16))
dev.off()

pdf("FigureA4-B.pdf",width=5,height=5)

ggplot(cz_taa, aes(estcertifiedwrks_n, petitioningworkers_n)) + xlab("Petitioning Workers, 1990-2007") + ylab("Certified Workers, 1990-2007") + geom_point() + geom_smooth() + theme(text = element_text(size=16))
dev.off()

# Appendix Figure A5 TAA Petitioners by Import Exposure per Worker in the 1990s (A) and 2000s (B) -- 2000-2010

pdf("FigureA5-A.pdf",width=7,height=5)

response_1990 <- mapPlot + geom_polygon(data=cz.df.1990, aes(x = long, y = lat, group = group, fill = cz.df.1990$response_qt9010), lwd = 1/10, colour="lavenderblush4") + scale_fill_continuous(name = "TAA Response to Import Exposure", low="white", high="gray30", guide="colorbar", limits=c(0.0,1.0), labels=c("0", "0.085", "0.313", "0.765", "770.04")) + theme(legend.text = element_text(size=10)) + coord_map(project="conic", lat0 = 30)

plot(response_1990)
dev.off()

pdf("FigureA5-B.pdf",width=7,height=5)

response_2000 <- mapPlot + geom_polygon(data=cz.df.2000, aes(x = long, y = lat, group = group, fill = cz.df.2000$response_qt9010), lwd = 1/10, colour="lavenderblush4") + scale_fill_continuous(name = "TAA Response to Import Exposure", low="white", high="gray30", guide="colorbar", limits=c(0.0,1.0), labels=c("0", "0.085", "0.313", "0.765", "770.04")) + theme(legend.text = element_text(size=10)) + coord_map(project="conic", lat0 = 30)

plot(response_2000)
dev.off()

# Appendix Figure A6 TAA-Certified Workers by Import Exposure in the 1990s (left) and 2000s (right)

pdf("FigureA6-A.pdf",width=7,height=5)

response_1990 <- mapPlot + geom_polygon(data=cz.df.1990, aes(x = long, y = lat, group = group, fill = cz.df.1990$cert_response_qt), lwd = 1/10, colour="lavenderblush4") + scale_fill_continuous(name = "Certified Workers per Import Exposure", low="white", high="gray30", guide="colorbar", limits=c(0.0,1.0), labels=c("0", "0.028", "0.174", "0.499", "701.12")) + theme(legend.text = element_text(size=10)) + coord_map(project="conic", lat0 = 30)

plot(response_1990)
dev.off()

pdf("FigureA6-B.pdf",width=7,height=5)

response_2000 <- mapPlot + geom_polygon(data=cz.df.2000, aes(x = long, y = lat, group = group, fill = cz.df.2000$cert_response_qt), lwd = 1/10, colour="lavenderblush4") +  scale_fill_continuous(name = "Certified Workers per Import Exposure", low="white", high="gray30", guide="colorbar", limits=c(0.0,1.0),  labels=c("0", "0.028", "0.174", "0.499", "701.12")) + theme(legend.text = element_text(size=10)) + coord_map(project="conic", lat0 = 30)

plot(response_2000)
dev.off()

# Appendix Figure A7 Change in TAA Responses to  Import Exposure, 1990-2007

pdf("FigureA7.pdf",width=7,height=5)

response_delta <- mapPlot + geom_polygon(data=cz.df.change, aes(x = long, y = lat, group = group, fill = cz.df.change$response_delta_qt), lwd = 1/10, colour="lavenderblush4") + scale_fill_gradient2(name ="Change in TAA Response", low="#B2182B",mid="white",high="#2166AC",midpoint=0.7603878, guide="colorbar",  limits=c(0.0,1.0), labels=c("-770.04", "-1.04", "-0.27", "0", "25.62")) + theme(legend.text = element_text(size=10)) + coord_map(project="conic", lat0 = 30)

plot(response_delta)
dev.off()

# Appendix Figure A8 Change in TAA Responses to  Import Exposure, 1990-2010

pdf("FigureA8.pdf",width=7,height=5)

response_delta <- mapPlot + geom_polygon(data=cz.df.change, aes(x = long, y = lat, group = group, fill = cz.df.change$response_delta_qt9010), lwd = 1/10, colour="lavenderblush4") + scale_fill_gradient2(name ="Change in TAA Response", low="#B2182B",mid="white",high="#2166AC",midpoint=0.7396122, guide="colorbar",  limits=c(0.0,1.0), labels=c("-770.04", "-1.03", "-0.25", "0", "85.15")) + theme(legend.text = element_text(size=10)) + coord_map(project="conic", lat0 = 30)

plot(response_delta)
dev.off()

# Appendix Figure A9 Import exposure in the 1990s (A),  2000s (B), and their change (C)

pdf("FigureA9-A.pdf",width=7,height=5)

exposure_1990 <- mapPlot + geom_polygon(data=cz.df.1990, aes(x = long, y = lat, group = group, fill = cz.df.1990$exposure_qt), lwd = 1/10, colour="lavenderblush4") + scale_fill_continuous(name = "Import Exposure per Worker ($1,000)", low="white", high="brown4", guide="colorbar", limits=c(0.0,1.0), labels=c("0", "0.11", "0.35", "0.93", "41.07")) + theme(legend.text = element_text(size=10)) + coord_map(project="conic", lat0 = 30)

plot(exposure_1990)
dev.off()

pdf("FigureA9-B.pdf",width=7,height=5)

exposure_2000 <- mapPlot + geom_polygon(data=cz.df.2000, aes(x = long, y = lat, group = group, fill = cz.df.2000$exposure_qt), lwd = 1/10, colour="lavenderblush4") + scale_fill_continuous(name = "Import Exposure per Worker ($1,000)", low="white", high="brown4", guide="colorbar", limits=c(0.0,1.0), labels=c("0", "0.11", "0.35", "0.93", "41.07")) + theme(legend.text = element_text(size=10)) + coord_map(project="conic", lat0 = 30)

plot(exposure_2000)
dev.off()

pdf("FigureA9-C.pdf",width=7,height=5)

exposure_delta <- mapPlot + geom_polygon(data=cz.df.change, aes(x = long, y = lat, group = group, fill = cz.df.change$exposure_delta_qt), lwd = 1/10, colour="lavenderblush4") + scale_fill_gradient2(name ="Change in Import Exposure", low="#2166AC",mid="white",high="#B2182B",midpoint=0.0719, limits=c(0,1.0), labels=c("-6.21", "0.21", "0.46", "1.01", "32.37")) + theme(legend.text = element_text(size=10)) + coord_map(project="conic", lat0 = 30)

plot(exposure_delta)
dev.off()

# Appendix Figure A10 TAA Petitioners in the 1990s (A), 2000s (B), and their change (C)

pdf("FigureA10-A.pdf",width=7,height=5)

taa_1990 <- mapPlot + geom_polygon(data=cz.df.1990, aes(x = long, y = lat, group = group, fill = cz.df.1990$taa_qt), lwd = 1/10, colour="lavenderblush4") + scale_fill_continuous(name = "TAA Petitioning Worker (%)", low="white", high="dodgerblue4", guide="colorbar", limits=c(0.0,1.0), labels=c("0", "0.026", "0.142", "0.366", "2.89")) + theme(legend.text = element_text(size=10)) + coord_map(project="conic", lat0 = 30)
plot(taa_1990)
dev.off()

pdf("FigureA10-B.pdf",width=7,height=5)

taa_2000 <- mapPlot + geom_polygon(data=cz.df.2000, aes(x = long, y = lat, group = group, fill = cz.df.2000$taa_qt), lwd = 1/10, colour="lavenderblush4") + scale_fill_continuous(name = "TAA Petitioning Worker (%)", low="white", high="dodgerblue4", guide="colorbar", limits=c(0.0,1.0), labels=c("0", "0.026", "0.142", "0.366", "2.89")) +  theme(legend.text = element_text(size=10)) + coord_map(project="conic", lat0 = 30)
plot(taa_2000)
dev.off()

pdf("FigureA10-C.pdf",width=7,height=5)

taa_delta <- mapPlot + geom_polygon(data=cz.df.change, aes(x = long, y = lat, group = group, fill = cz.df.change$taa_delta_qt), lwd = 1/10, colour="lavenderblush4") +  scale_fill_gradient2(name ="Change in TAA Applications", low="#B2182B",mid="white",high="#2166AC",midpoint=0.44182825, guide="colorbar",  limits=c(0.0,1.0), labels=c("-2.31", "-0.07", "0", "0.14", "1.94")) +  theme(legend.text = element_text(size=10))
taa_delta <- taa_delta + coord_map(project="conic", lat0 = 30)

plot(taa_delta)
dev.off()